#include<iostream>
using namespace std;
#include<cmath>

/*
 * Nevill algorithm
 *
 * Given n data pointes, in (x_data, y_data), and an arbitrary x value, the function Nevill outputs y=P(x), and error estimation dy
 *
 * If the datapoints are too close, round off error can be significant, therefore the function refuses the task
 */
void Nevill(int n, double *x_data, double *y_data, double x, double *y, double *dy){
	//check if data points are too close
	for(int i=0;i<n;i++){
		for(int j=i+1;j<n;j++){
			if( fabs(x_data[i]-x_data[j]) < 1E-6 ){
				cout<<" datapoints' x values are too close.\n";
				exit(1);
			}
		}
	}

	double * tempP = new double [n];
	for(int i=0;i<n;i++) tempP[i] = y_data[i];// initially P_0, P_1, \cdots, P_{n-1}
	//use the recurrence relations to derive the final polynomial
	for(int i=1;i<n;i++){
		for(int j=i;j<n;j++){
			tempP[j] = (x-x_data[j])*tempP[j-1] + (x_data[j-i] - x)*tempP[j];
			tempP[j] /= x_data[j-i] - x_data[j];
		}
	}

	*y = tempP[n-1];
	*dy = 0;
	if(n>1) *dy = fabs(tempP[n-1] - tempP[n-2]);
	delete [] tempP;
}

int main(){

	int n;
	double x,y,dy;

	double *x_data = new double [n];
	double *y_data = new double [n];

	FILE *fp=fopen("Nevill_input.txt","r");
	if(fp==NULL){
		cout<<"error: failed to open Nevill_input.txt.\n";
		exit(1);
	}

	n=0;
	while( fscanf(fp, "%*c%*[^\n]") != EOF ){
		n++;
	}
	cout<<"n="<<n<<endl;
	
	fp=fopen("Nevill_input.txt","r");
	for(int i=0;i<n;i++){
		fscanf(fp, "%lf %lf", &x_data[i], &y_data[i]);
	}
	fclose(fp);

	double xmin = x_data[0], xmax = x_data[0];
	for(int i=0;i<n;i++){
		if( xmin > x_data[i] ) xmin = x_data[i];
		if( xmax < x_data[i] ) xmax = x_data[i];
	}
	xmin -= 2;
	xmax += 2;

	int num_grid = 100;
	double step = (xmax - xmin)/num_grid;

	fp=fopen("Nevill_curve.txt","w");
	if(fp==NULL){
		cout<<"error: failed to open Nevill_curve.txt.\n";
		exit(1);
	}
	for(int i=0;i<num_grid;i++){
		x = xmin + step*i;
		Nevill(n,x_data,y_data,x,&y,&dy);
		fprintf(fp,"%lf   %lf   %lf \n", x, y, dy);
	}
	fclose(fp);

	return 0;
}
